rm(list = ls())
library(tidyverse)
library(ggplot2)
# parameters set
A <- 1
dlt <- 0.1
sgm <- 0.2
theta <- 0.36
n <- 0.02
eps <- rnorm(120, sd = 0.2)

# 迭代. 公式(1.2)
compk <- function(A, dlt, sgm, theta, n, eps){
  k <- rep(2.2,120) # 初值2.2
  for (i in 2:120) {
    k[i] <- ((1-dlt)*k[i-1] + sgm*A*exp(eps[i-1])*k[i-1]^theta)/(1+n)
  }
  return(data.frame(k = k))
}

# 模拟三次
set.seed(1024)
eps <- lapply(rep(0.2,3),rnorm, n = 120, mean = 0)
picdata <- lapply(eps, compk, A = A, dlt = dlt, sgm = sgm, theta = theta, n = n) %>% bind_cols()
names(picdata) <- c('k1','k2','k3')
picdata$time <- 1:120

# 画图
ggplot(picdata, aes(x = time, y = k1)) + geom_line(linetype = 1) +
  geom_line(aes(y = k2), linetype = 2) +
  geom_line(aes(y = k3), linetype = 9) + labs(y = 'k') +
  theme_bw()
ggsave('../EconomicGrowthLecture/solow1.png')

# 对数线性化的迭代. 公式(1.4)
B <- (1+theta*n-dlt*(1-theta))/(1+n)
C <- (dlt+n)/(1+n)

compk_ln <- function(B, C, eps){
  k <- rep(0,120) # 因为是对数差分的情况，所以初值为0
  for (i in 2:120) {
    k[i] <- B*k[i-1] + C*eps[i-1]
  }
  return(data.frame(k = k))
}


# 画图
picdata$k_ln <- compk_ln(eps = eps[[1]],B = B, C = C) %>%
  unlist()
picdata$k_ln <- 2.2*exp(picdata$k_ln) # 由对数差分返回原值

ggplot(picdata, aes(x = time, y = k1)) + geom_line(linetype = 1) +
  geom_line(aes(y = k_ln), linetype = 2) + labs( y = 'k') +
  theme_bw()
ggsave('../solow2.png')
